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1. Introduction 

The free Dirac equation predicts that an elementary fermion has an intrinsic magnetic moment 
proportional to its spin S, 




(1.1) 



where the Lande g-factor is exactly 2, e is the fundamental electric charge, and m the fermion's mass 
(H = c = 1). In the presence of interactions, g receives radiative corrections (Fig. 1), and these can 
be computed order-by-order in weak coupling perturbation theory in the electromagnetic coupling 
a = e 2 /An. As discussed below, in addition, the hadronic corrections require non-perturbative 
methods or experimental data. 

The fundamental vector interaction of the fermion with an external electromagnetic field is 
given by the matrix element of the electromagnetic current between incoming and outgoing fermion 
states of momentum and spin {p,s} and {p' ,s'}, 

(p',s'\Y^Y\P,s) = u(p',s') ^(Q^ + i^^-ir^^Qv^) u(p,s), (1.2) 

where u(p,s), il(p',s') are spinors, and Q v = (p — p') v is the (space-like) momentum transferred 
to the fermion by the photon. The form of Eq. (1.2) is dictated by Lorentz-invariance, P and C 
symmetries, and the Ward-Takahashi identity; hence the form factors f\(<2 2 ) and /^(G 2 ) contain 
all information on the fermion's interaction with the photon. In particular, in the limit Q 2 — > 0, 
F\ (0) is the charge of the fermion in units of e, and 

F 2 (0) = ^ (1.3) 

is the anomalous magnetic moment. In the following 7*2(0) = is referred to as the muon 
"anomaly". 

The muon anomaly provides one of the most stringent tests of the Standard Model because it 
has been measured to fantastic accuracy (0.54ppm) [1] and calculated to even better precision [2, 
3, 4]. The difference between the two is reported to range between 249(87) x 10~ n and 287(80) x 
10~ n , or about 2.9 to 3.6 standard deviations [2, 3, 4]. The Standard Model precision is attained by 
calculating contributions (in orders of a) from QED (a 5 ) [2], electroweak (EW) (a 2 ) [5], and QCD 
(a 3 ). The latter represent the largest uncertainty. The lowest order (a 2 ) contribution arises from 
the hadronic vacuum polarization (HVP), Yl(Q 2 ), while the next-to-leading contribution is from 
hadronic light-by-light (HLbL) scattering and the HVP accompanied by an additional photon. The 
various contributions to are summarized in Tab. 1 and depicted in Fig. 1. It is well known 
that the discrepancy between theory and experiment is large in magnitude, roughly two times the 
EW contribution, and has important consequences for new physics. Because of its central role in 
constraining beyond the Standard Model (BSM) theories, the muon anomaly will be measured to 
even better accuracy in new experiments at Fermilab (E989) and J-PARC (E34). The goal of E989 
is to reduce the total error to 0. 14 ppm. A brief summary of the importance of (g — 2)^ to BSM 
physics can be found in [7]. 
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Table 1: Standard Model contributions to the muon anomaly. The QED contribution is through a 5 , EW 
a 2 , and QCD a 3 . The two QED values correspond to different values of a, and QCD to lowest order (LO) 
contributions from the hadronic vacuum polarization (HVP) using e + e~ —> hadrons and T — > hadrons, higher 
order (HO) from HVP and an additional photon, and hadronic light-by-light (HLbL) scattering. 



QED 




11658471.8845 (9)(19)(7)(30) x 1(T 10 


[2] 






11658471.8951 (9)(19)(7)(77) x 1(T 10 


[2] 


EW 




15.4(2) x 1(T 10 


[5] 


QCD 


LO (e+e-) 


692.3 (4.2) x 1(T 10 , 694.91 (3.72) (2.10) x 1(T 10 


[3,4] 




LO (T) 


701.5(4.7) x 10~ 10 


[3] 




HO HVP 


-9.79(9) x 10~ 10 


[6] 




HLbL 


10.5(2.6) x 10- 10 


[9] 



The HVP contribution to the muon anomaly has been computed using the experimentally 
measured cross-section for the reaction e + e~ — > hadrons and a dispersion relation to relate the real 
and imaginary parts of Yl(Q 2 ). The current quoted precision on such calculations is a bit more than 
one-half of one percent [3, 4]. The HVP contributions can also be calculated from first principles 
in lattice QCD [8]. While the current precision is significantly higher for the dispersive method, 
lattice calculations are poised to reduce errors significantly in next one or two years. These will 
provide important checks of the dispersive method before the new Fermilab experiment. Unlike 
the case for (HVP) , (HLbL) can not be computed from experimental data and a dispersion 
relation (there are many off-shell form factors that enter which can not be measured). While model 
calculations exist (see [9] for a summary), they are not systematically improvable. A determination 
using lattice QCD where all errors are controlled is therefore desirable. 

In Sec. 2 we review the status of lattice calculations of ^(HVP). Section 3 is a presentation 
of our results for (HLbL) computed in the framework of lattice QCD+QED. Section 4 gives our 
conclusions and outlook for future calculations. 




Figure 1: Representative diagrams, up to order a 3 , in the Standard Model that contribute to the muon 
anomaly. The rows, from to top to bottom, correspond to QED, EW, and QCD. Horizontal solid lines 
represent the muon, wiggly lines denote photons unless otherwise labeled, other solid lines are leptons, 
filled loops denote quarks (hadrons), and the dashed line represents the higgs boson. 
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2. Hadronic Vacuum Polarization 

It has long been known that (HVP) can be computed from e + e~ — > hadrons [10, 1 1]. Later it 
was realized that % — > hadrons could also be used to improve the result in a limited, but important, 
range of Q 2 [12]. For example, the most recent studies find [4] ^(HVP) = (694.91 ± 3.72 ± 
2.10) x 10~ 10 (e+O while [3] quotes a M (HVP) = (692±4.2) x 10~ 10 (e + e~) and (701.5 ±4.7) x 
10~ 10 (including t). These determinations are the targets for lattice calculations. Note the result 
including X data is about 2 standard deviations larger than the pure e + e~ contribution, and reduces 
the discrepancy with the Standard Model to 2.4 standard deviations [3]. The former requires isospin 
corrections which may not be under control. Alternatively, p — / mixing may explain the difference 
and bring the T-based result in line with that of e + e~ [13]. The lattice QCD calculations, as we 
discuss below, are not yet at the level of precision to help resolve this issue. 

In [8] it was shown that fl^(HVP) can be calculated instead using lattice QCD, 

(Qf\2 f°° 
-) I dQ 2 f(Q 2 )tl(Q 2 ). (2.1) 

(a similar formula was found much earlier in the context of large N QCD [14]). Here Q 2 is the 
Euclidean momentum-squared in the one-loop QED integral (Fig. 2) that runs through the renor- 
malized (once-subtracted) vacuum polarization, tl(Q 2 ), and f(Q 2 ) is a known analytic function 
that diverges in the limit Q 2 — > 0. Hence the integrand is dominated by the low momentum region 
around m 2 . The HVP is non-perturbative in this region, so the lattice provides a natural framework 
for the calculation of ^(HVP). The importance of the low momentum region (and light quark 
masses) implies large volumes are necessary. Since lattice results are determined at discrete values 
of Q 2 , a smooth interpolation is necessary in order to utilize Eq. (2.1). tl{Q 2 ) is also cut-off at 
\Q\ ~ I /a, so the high-momentum region can be handled with continuum perturbation theory [15], 
though this contribution is small [8], on the order of a percent for present lattice calculations. 





„. - „. , j i j ■ i Figure 3: Lowest order hadronic contribution to 

tigure 2: The lowest order hadronic contnbu- 

, . „, , , , the muon anomaly from two quark loops con- 

tion to the muon anomaly. The blob represents 

.. . . ,. . . nected by gluons. The contribution is flavor- 

all possible intermediate hadronic states. 

symmetry and Zweig suppressed. 

First we briefly review the published lattice results for (HVP) which are summarized in 
Tab. 2 (and have not changed since last year's meeting [16]). We then update them from this year's 
meeting and discuss new developments. We discuss only results using dynamical fermions as it is 
known that the quenched results significantly underestimate the dispersive ones [8, 17]. The quoted 
eiTors on all of the results are still somewhat crude, and nowhere complete. The first dynamical 
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results were obtained by Aubin and Blum using 2+1 flavors of improved staggered quarks [18]. The 
statistical errors are 2-3%, but systematic errors are not controlled: there is only one lattice spacing 
so far, and the fits to Yl(Q 2 ), which are based on a resonance model+chiral perturbation theory for 
staggered quarks, are very sensitive to the low-momentum region. In addition, the vector mass in 
the pole is fixed to a value determined in a separate analysis. While a strength of the calculation 
is that Goldstone pion masses as light as m % « 220 MeV are used, the simple linear and quadratic 
extrapolations in the quark mass are problematic since the p -meson is stable in the calculation (as 
is true in all current lattice calculations, though at least one calculation where the p is unstable is 
underway [19]). 

The UKQCD group has computed ^(HVP) using 2+1 flavors of domain wall fermions with 
several lattice spacings, volumes, and quark masses (down to m n 170 MeV) [20]. They also 
developed a momentum fit based on two vector-particle poles in n(<2 2 ). In order to stabilize the 
fit, the lowest pole is fixed to the vector mass determined on the same lattices but from a separate 
analysis. The lattice spacing dependence can not be resolved outside relatively large statistical 
errors, so all results are fit together in a lattice-spacing-independent fit. An alternative quark mass 
fit [21] provides an estimate of the chiral extrapolation error which is about 5%, the same size as 
the statistical error. 

The Mainz group uses two flavors of Wilson fermions with several volumes, relatively heavy 
pion masses (down to 280 MeV), a quenched strange quark, and several lattice spacings, but quotes 
results for a^(HVP) from only one lattice spacing [22]. A partial systematic error of about 10% 
coming from cut-off effects is estimated by re-analyzing the data with the lattice spacing shifted by 
one statistical standard deviation. The alternative chiral exptrapolation in [21] is also tried here, and 
no significant difference is found. The statistical error is also about 10%. An innovation first used 
here that will likely be used in future calculations is twisted-boundary conditions for the quarks, so 
lower values of the momenta than allowed by periodic boundary conditions are accessible. 

Finally, the ETM collaboration has computed (HVP) in two flavor QCD with twisted-mass 
Wilson fermions [21]. They employ several lattice spacings, relatively small volumes, and heavy 
pion masses. The alternative chiral extrapolation mentioned above was proposed here and used to 
quote central values: The kernel in Eq. (2.1) is rescaled by a reference hadronic observable, e.g., 
my or fy. Nominally, this is designed to lead to smaller quark mass dependence, but in this case 
simply yields an overall shift in the data towards larger values of (HVP) since the quark masses 
used are far from the chiral regime. As in the other studies, effects of volume and lattice spacing 
are not detectable within relatively large statistical errors, so all data is included in a lattice-spacing 
and volume independent fit. The final quoted error is simply the statistical uncertainty from this fit 
and is about 3%. 

At this year's conference, two of the above groups presented preliminary results of improved 
calculations and methods. Jager presented results for the Mainz group with updated statistics and 
lighter quark masses (m n > 195 MeV) [23]. Hotzel presented preliminary 2+1+1 results from 
the ETM collaboration [24] in which the lightest pion mass is down to about 270 MeV. From the 
dispersive calculations the contribution from the charm quark is estimated to be about the same as 
the HLbL contribution, a^(HVP, charm) « tty(HLbL) « 10 X 10~ 1() , and they find about double 
that amount. 

Golterman presented a new model-independent method for fitting the Q 2 dependence of H(Q 2 ) 
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Table 2: Published lattice QCD results for a^{HWP)(x 10 10 ). 





N f 


errors 


action 


group 


713(15) 


2+1 


stat. 


Asqtad 


[18] 


748(21) 


2+1 


stat. 


Asqtad 


[18] 


641(33)(32) 


2+1 


stat., sys. 


DWF 


[20] 


572(16) 


2 


stat. 


TM 


[21] 


618(64) 


2+1 1 


stat., sys. 


Wilson 


[22] 



based on Pade approximants [25, 26]. The central idea is that since fl(Q 2 ) is given by a positive 
spectral function, it is a so-called Steiljes function, obeying certain constraints, and can be obtained 
from a converging series of Pade approximates evaluated on a sequence of Q 2 points. The Pade 
approximates furnish a natural, model-independent, fit function for H(Q 2 ). Unfortunately, it is 
found that the new method disagrees outside of statistical errors with the vector-meson-dominance 
(VMD) fit function used earlier in [18]. This discrepancy can be traced to the aforementioned 
sensitivity of (HVP) to the low momentum region. In Fig. 4 we show the two fits which lead to 
(HVP) = 350(8) and 413(8), respectively. Both fits are good in the sense that they go through the 
data points and have acceptable % 2 values (the VMD fit is uncorrelated, i.e. the covariance matrix 
is diagonal), yet they give quite different results. The statistical errors on H(Q 2 ) are not small 
enough to differentiate the two fits, leading to an uncomfortably large systematic error. To make 
progress they must be reduced. In general, using higher Q 2 points where the errors are smaller 
leads to systematically low fits in the low momentum region [18]. Relatedly, uncorrelated fits or 
fits with more free parameters mitigate the effect, but at the cost of larger statistical errors. 

One way to reduce the statistical errors on Tl(Q 2 ) in the critical low momentum region is to use 
twisted boundary conditions to fill in this region [22, 23] . This is not enough however, because the 
errors on each point are still too large. A new method described at the meeting by Shintani [27, 28] 
that addresses this problem goes under the name "All Mode Averaging" (AMA). In Fig. 5 we show 
the HVP computed using the AMA technique. Computer resources required to obtain the same 
statistical error on Yl{Q 2 ) in the range < Q 2 < 1.0 GeV 2 were reduced by factors of 2.6-20 over 
the conventional method. AMA should be even more effective as the quark mass is reduced and 
the volume is increased. 

So far we have only discussed the contribution to (HVP) arising from the diagram shown 
in Fig. 2, i.e., from a single quark loop. In 2+1 flavor QCD and the SU (3) flavor-symmetric limit, 
this is all there is. However flavor or isospin symmetry breaking gives rise to contributions from 
diagrams with two quark loops connected by gluons, as shown in Fig. 3. While this contribution 
is flavor-symmetry and Zweig suppressed, it is expected to contribute on the 1% level, possibly 
more, so must eventually be included. In [29], using chiral perturbation theory, the two-quark-loop 
("disconnected") diagram was shown to account for 10% of the connected contribution. However, 
it is known that the pion contribution is only a small part of the total [18]. The two-quark-loop 
contribution was calculated in [21], but with a statistical precision that was not significant. 
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Figure 4: The hadronic vacuum polarization 
computed on a 2+1 flavor Asqtad ensemble from 
the MILC collaboration [25, 26]. The solid 
line denotes a fit to the [1,1] Pade approximant 
while the dashed line is from a vector-meson- 
dominance fit. 



0.2 
0.18 
0.16 
0.14 
0.12 

0.1 
0.08 
0.06 




4 6 
Q 2 (GeV 2 ) 



10 



Figure 5: The hadronic vacuum polarization for 
a single flavor computed on a 2+1 flavor Asqtad 
ensemble from the MILC collaboration (m K = 
315 MeV) using AMA (circles) [28]. Original 
HVP calculation (squares) from [30]. 



3. Light by Light Amplitude 

The hadronic light-by-light (HLbL) scattering amplitude is shown in Fig. 6. Unlike the HVP 
contribution, there is no dispersive method of calculation for (HLbL) . Model estimates put the 
value at ^(HLbL) = 10.5(2.6) x 10~ 10 [9]. The uncertainty quoted for the estimate is part of the 
"Glasgow consensus" and is not attributable to any one calculation. The errors on the model cal- 
culations can not be systematically reduced, and therefore it is desirable to calculate the amplitude 
directly, using lattice QCD. Notice from Tab. 1 that the errors on (HVP) are roughly two times 
as large as the quoted model error. Some uncertainties in the literature are larger, for example, 
see [31]. In any case, as the error on the HVP contribution is reduced, the HLbL error will come to 
dominate the theory error. 

The diagram in Fig. 6 can, in principle, be computed from the correlation of four electromag- 
netic currents in QCD. Since the aim is to compute the corresponding two-loop QED integral, and 
one of the electromagnetic currents corresponds to the external photon whose momentum should 
be taken to zero, two independent momenta are needed for the three remaining vertices due to 
momentum conservation. This is naively a (N t x N^) 2 difficult problem. A simpler, but useful 
exercise, would be to compute the four-point correlation function at some number of fiducial mo- 
menta points to check the model calculations. One can also compute a part of the amplitude, like 
71 — > y* 7* (where at least on photon is off-shell) [32] which is also useful for model calculations. 

A different approach is to compute the entire amplitude, including the muon and photons, 
on the lattice [33]. One can either treat the photons in a non-perturbative framework, the proce- 
dure that will be followed here, or by treating QED perturbatively but still within the lattice QCD 
framework. The basic method has been laid out in [33] and is shown schematically in Fig. 8. The 
idea is to compute correlation functions in a combined photon plus gluon field configuration as 
first done in [34]. Start by forming a quark loop with two electromagnetic vertices, one for the 
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Figure 6: The hadronic light-by-light scattering 
amplitude. The quark loop represents all pos- 
sible ("connected") intermediate hadronic states 
connected to the horizontal muon line by three 
photons. The fourth photon represents the exter- 
nal electromagnetic field. 




Figure 7: A disconnected quark loop diagram 
that is part of the hadronic light-by-light scatter- 
ing amplitude. The quark loops are connected by 
gluons. 




Figure 8: Non-perturbative method for calculating the light-by-light amplitude. The full correlation function 
with one internal photon inserted by hand is computed in QCD+QED (left). The subtraction term (right) is 
constructed from separate ensemble averages in QCD+QED and QED. The loop represents the quark (or 
lepton in pure QED) and the horizontal line is for the muon. The difference of the two diagrams yields the 
light-by-light scattering amplitude. 



external photon and another for one of the internal ones. The loop is connected to the muon line 
by a free lattice photon propagator, so the correlation function is constructed from a two-point 
quark loop and a three point muon line, where the incoming and outgoing muons have momenta 
p and p', respectively, and is averaged over the combined QCD+QED gauge field. As it stands, 
the correlation function has a leading contribution proportional to a (which holds even if the QED 
configurations are omitted). To extract the HLbL amplitude, which is 0(a 3 ), a non-perturbative 
subtraction is necessary. Typical higher order terms in the correlation function are shown in Fig. 9. 
The subtraction term is constructed in exactly the same way, except at the last step, the loop and 
line are averaged separately over QCD+QED and QED configurations, respectively. The subtrac- 
tion has exactly the same terms as the original correlation function, through 0(a 3 ), except the 
desired HLbL term which is absent since only one photon connects the quark loop and the muon 
line (see Fig. 8). The difference of the two yields the HLbL amplitude. The key point in practice 
is that both the correlation function and the subtraction term use exactly the same QED and QCD 
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configurations, so they are highly correlated. By using suitable projectors on the subtracted cor- 
relation function, linear combinations of the form factors in Eq. (1.2) are obtained (i.e., the Sachs 
form factors Ge and Gm), from which F2 is found. 



Figure 9: Leading and higher order diagrams in the expansion of the correlation and subtraction functions 
shown if Fig. 8. The diagram on the left is (a), the next two are a 2 , and the last two are a 3 . 

In practice, the point-split vector current, which is conserved at non-zero lattice spacing, is 
used at the internal vertices for the muon line and quark loop so that Ward-Takahashi identities are 
satisfied exactly, up to numerical precision, on each configuration. These vertices are each Fourier 
transformed to 4-momentum space, and both the full correlation function and the subtraction cor- 
relation function are multiplied by the photon propagator and summed over this 4-momentum. For 
simplicity, a local current (YYnY) which must be renormalized is used for the external vertex. 

3.1 QED contribution 

In order to check the method we started with the calculation in pure QED where the answer 
is well known in perturbation theory [35, 36]. The calculation was done in quenched non-compact 
QED, in the Feynman gauge, using domain wall fermions (DWF) [37]. The lattice size was 16 3 x 32 
with L s = 16 sites in the extra 5th dimension. The muon mass in this test case was relatively 
large, 0.4 in lattice units, and to enhance the signal the electric charge was set to e = 1.0 which 
corresponds to a = \/4n instead of 1/137. Here and in the following we always use kinematics 
where the incoming muon is at rest. 

The first attempt was made with the same heavy mass in the lepton loop as for the muon line, 
which amounts to computing the electron g — 2. Using a single point source at the external vertex 
and relying on plane wave sources for the incoming and outgoing muons to conserve momentum, 
we were unable to attain a signal (Fig. 10), though the statistical error, obtained from a few hundred 
configurations, was the size of the signal expected from perturbation theory. Taking advantage of 
the logarithmic enhancement due to the lighter electron mass, ln(m^/m e ) [36], the loop mass was 
reduced to 0.01 . Accordingly, the signal should have increased by an order of magnitude. The result 
is shown in the right panel of Fig. 10; a clear signal is observed. The form factor F2 was computed 
only at the lowest non-trivial momenta, 2jt/16, and was not extrapolated to zero. The average of 




F 2 (Q 2 ), 



(3.1) 



(3.2) 




9 



Hadronic corrections to the muon g — 2 from lattice QCD 



T. Blum 



the three points in the "plateau" shown in Fig. 10 give F 2 = 3.96(70) x 10 -4 = 24.4(4.3) (a/n) 3 
while perturbation theory gives about 10(a/7l) 3 for ^(0). The calculation was repeated on a 
24 3 x 32 lattice to investigate finite volume effects. Again, for the lowest non-trivial momentum, 
we find F 2 = 1.19(32) x 10~ 4 = 7.32(1.97)(a/7r) 3 , roughly consistent with perturbation theory. 
The renormalization factor of the local vector current, inserted at the external vertex, was not 
included, but its effect which is 0(a) should be small compared to other uncertainties. 



• 500 configs (Coulomb) 

* 500 configs (Coulomb) 
+ 300 configs (Feynmann) 



- b 



0001 
0.0009 
0.0008 
0.0007 
0.0006 
' 0.0005 
0.0004 
0.0003 
0.0002 
0.0001 




Figure 10: The form factor Fi(Q 2 ) for light-by-light scattering in QED [37], evaluated at the lowest non- 
nivial lattice momenta, q = 2n/\6. The mass of the lepton in the loop is 0.4 (left) and 0.01 (right). t op is the 
time slice of the external electromagnetic vertex. 



3.2 QCD contribution 

The inclusion of QCD into the light-by-light amplitude is straightforward: simply multiply 
the £7(1) gauge links with SU (3) links to create a combined photon and gluon configuration [34], 
and follow exactly the same steps, using the same code, as described in the previous sub-section. 
We use one QED configuration per QCD configuration, though different numbers of each could 
be beneficial and should be explored. Since the photons are quenched in this first calculation, the 
valence quarks are charged but not the sea. We discuss the implications of this below. 

For QCD we start with the 16 3 x 32 x 16, 2+1 flavor, DWF+Iwasaki ensemble generated 
by the RBC/UKQCD collaboration [38]. The light sea quark mass is 0.01 , the strange 0.032, 
and the residual mass m res = 0.00308(4). The gauge coupling is jS = 2.13 and inverse lattice 
spacing aT x = 1.73(3) GeV. The corresponding pion mass is m% = 422 MeV, and the muon mass 
is = 0.4 2 is the same as in the pure QED case. The external electromagnetic vertex is inserted 
on time slice t op = 6 and the incoming and outing muons are created and destroyed at t = and 
12, respectively. The incoming muon is again at rest while the outgoing muon has one unit of 
momentum (we average results over all three directions, and ±p). Finally, e = 1 as before. While 
these unphysical parameters may induce relatively large systematic errors, the main aim of this 
study is to show that the statistical errors can be controlled and that the method works. After 
demonstrating this, the systematic errors can be addressed. 

In the QED case, setting m e <C lead to a large enhancement of the signal to noise ratio. The 
same effect is not expected here since m n ~ m^. Even in the physical case, m K /m^ 1.33 which 

2 Using the QCD lattice spacing, this corresponds to 692 MeV compared to the physical mass, = 105.658367(4). 
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may partly explain why the HLbL amplitude is about 200 times smaller than the QED amplitude. 
To reduce the statistical error we use a point source at the external vertex of the quark loop every 
4 spatial sites in each direction, for a total of 64 propagator calculations per configuration (100 
combined-configurations were used). 

Following the same procedure as in the pure QED case, we find F2(Q 2 = 0.38 GeV 2 ) = 
(—15.7 ±2.3) x 10~ 5 . While the signal is robust, the magnitude is between 5-10 times larger 
than the model estimates for a^(HLbL), and even has the opposite sign. There are several reasons 
why this might be the case. First, models are not expected to be valid in the range of masses used 
here. Since the pion mass is heavy, the leading pion pole term is probably not dominant and the 
sub-leading terms can be important [9]. Some of them have opposite signs [9]. Further, the sub- 
leading terms are known to be proportional to m 2 [9], enhancing the effect of the large muon mass 
used here. 

An important check that the subtraction is working is to vary e and see that the amplitude 
changes in the expected proportion compared to the case where e = 1. The same non-compact 
QED configurations are used in each case; e is varied only when constructing the exponentiated 
gauge-link, U^(x) = expieA^(x). Thus the ratio of form factors, and hence a dependence, can be 
determined very accurately. Since one photon is inserted explicitly, and the charge at this vertex 
is not included in the lattice calculation, the subtracted amplitude should behave like ~ e 4 . Using 
e = 0.84 and 1.19, the changes relative to e = 1 should be 0.5 and 2.0, respectively. This is precisely 
what is observed. In addition, there is no detectable change in the unsubtracted amplitudes. 

The next step, which is not too costly, is to reduce only the muon mass. With =0.1 and 
everything else unchanged, except that the number of configurations was increased to 306, we 
find F 2 {Q 2 = 0.19GeV 2 ) = -2.21 ±0.75 x 10~ 5 . This reduction is consistent with the already 
mentioned leading m 2 dependence. Note the value of Q 2 is also a factor of two smaller because of 
the smaller muon mass. 

Our next step was to increase the lattice size from 16 3 x 32 to 24 3 x 64 and reduce the quark 
mass to 0.005 (also an RBC/UKQCD 2+1 flavor, DWF+Iwasaki ensemble [39, 40]). The pion mass 
measured on this ensemble is 329 MeV, and volume is increased by more than a factor of three. 
The increased lattice size and smaller quark mass lead us to use the all mode averaging (AMA) 
technique [28] to maintain large statistics at an affordable cost. Besides the exact calculation which 
was done using a single point source on 20 configurations, the approximation was computed using 
400 lowmodes of the even-odd preconditioned Dirac operator and 216 point sources computed 
with stopping residual 10 on 157 configurations 3 . The external vertex was inserted on time 
slice t op = 5 with incoming and outgoing muons at and 9, respectively. The incoming muon is 
at rest, while the outgoing has three-momenta in units of 2n/L of (±1,0,0) and (±1 , ±1,0) and 
permutations. We also include the vector current renormalization (in pure QCD) from [40]. 

In Fig. 1 1 we show the stability of the signal for /^(G 2 ) for the two lowest non- trivial momenta 
(Q 2 = 0.11 and 0.18 GeV 2 ) as more configurations are added to the ensemble average (20 to 157, 
roughly doubling at each step). The statistical errors scale like l/VN, where N is the number of 
measurements. The central value is stable, and has magnitude comparable in the size to model 

3 This is about double the statistics as were presented at the meeting. 
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estimates (and statistical errors only 2-3 times larger). The sign is also the same. 

F 2 (0.18GeV 2 ) = (0. 1 15 ± 0.044) x 
F 2 (0.11GeV 2 ) = (0.014 ± 0.063) x 
a M (HLbL/model) = (0.084 ± 0.020) x 

Several comments are in order. The signal for the larger Q 2 point is about 3 standard devia- 
tions while the lowest Q 2 point is zero within the statistical error. The former has more momenta 
combinations to average over, 12 versus 6, which is reflected in the errors. The fact that the sign is 
consistent with model estimates (and is the same as for pure QED) may be due to the reduced pion 
mass (pion pole term becoming dominant), or the larger volume, or a combination of the two. In 
the AMA procedure [28] the expectation value of an operator is given by 

{&) = Kest> + ^£(^approx,g>, (3.6) 

where No is the number of measurements of the approximate observable, and "rest" refers to the 
contribution of the exact observable minus the approximation, evaluated for the same conditions. In 
the present calculation, the central value and errors appearing in Eqs. (3.3) and (3.4) are completely 
dominated by the second term on the r.h.s. of Eq. (3.6). 

We briefly comment on the absence of sea-quark electric charge in our calculation. The QED 
quenching of the sea-quarks is equivalent to omitting diagrams like the one shown in Fig. 7, and 
variants with three quark loops, or ones generated by interchanging loops containing the external 
vertex. In principle these could be calculated in the quenched setup, by introducing one and two 
additional valence quark loops. However it is straightforward to include them by repeating the 
calculation above, but on dynamical QCD+QED configurations, or equivalently, by re-weighting 
the quenched QED configurations [41]. We also note that there is no reason, a-priori, to expect 
that the size of the omitted diagrams is any smaller than what has already been calculated, except 
possibly by Zweig suppression. 

4. Conclusions 

The muon anomalous magnetic moment continues to provide one of the most precise tests 
of the Standard Model of particle physics and often places important constraints on new theories 
beyond the Standard Model. With new experiments planned at Fermilab (E989) and J-PARC (E34) 
that aim to improve on the current 0.54 ppm measurement at BNL by at least a factor of four, it 
will continue to play a central role in particle physics for the foreseeable future. 

Owing to the non-perturbative nature of QCD, the hadronic corrections to the muon g — 2 are 
the largest source of error in the Standard Model calculation. These errors must be reduced to 
leverage the new experiments [7]. The hadronic corrections enter at order a 2 through the hadronic 
vacuum polarization (HVP) and a 3 through hadronic light-by-light (HLbL) scattering as well as 
higher order HVP contributions. Both are being calculated using lattice QCD. 
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Figure 11: F2(Q 2 ) computed on the 24 3 , m n = 329 MeV ensemble. Lowest two momenta only. The 
statistical errors scale with 1/VN where N is the number of measurements which range from 20 to 157, 
roughly doubled at each step. The central value is also stable, and the size of model estimates. 



While the HVP contribution to the muon anomaly has been precisely computed to the 0.6% 
level using experimental measurements of e + e~ — > hadrons and z — > hadrons, lattice QCD calcula- 
tions serve as an important independent check of these results. At the moment, statistical errors on 
lattice calculations of a^(HVP) are at about the 3-5% level, but important systematic errors remain. 
Most significant is that for light quark masses the errors on the low-momentum region of Yl{Q 2 ) 
are not small enough, nor are there sufficient points available in the crucial region, Q 2 ~ m 2 ^, to 
adequately estimate ^(HVP). Quark masses are too heavy (or errors are still too large for light 
masses) and fits are model-dependent. The good news is that all of these points are being addressed 
in the newest calculations. Lattice calculations using model independent fit functions, noise reduc- 
tion techniques, twisted boundary conditions, charm quarks, and physical light quark masses on 
large lattices are underway: Big error reductions over the next one to two years are not only possi- 
ble, but likely. Finally, to get to the 1 % level, or better, disconnected diagrams and isospin breaking 
effects must be incorporated to complete the calculation. 

The HLbL contribution to the muon anomaly can not be computed using data from experiment 
and a dispersion relation as for the HVP contributions. It is expected to dominate the theory error 
as the HVP error is reduced (either from more experimental data, lattice calculations, or both). We 
have presented encouraging first results for the part of the HLbL contribution where there is a single 
quark loop in the corresponding diagram. Calculations in pure QED were used to successfully 
check the method. The calculation is numerically intensive, but logically straightforward. Much 
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effort is still needed to reduce statistical errors, extrapolate to zero momentum transfer, and many 
systematic errors remain uncontrolled. The calculation is quenched with respect to QED, so the sea 
quarks are not charged, and potentially large contributions are missing. This can be fixed by using 
dynamical QCD+QED gauge configurations, or by re-weighting the quenched QED ensembles. 
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